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ABSTRACT 

Polynomial interpolation is an essential subject in numerical analysis. Dealing with a 
real interval, it is well-known that even if f(x ) is an analytic function, interpolating at 
equally spaced points can diverge [Davi75]. On the otherhand, interpolating at the zeroes 
of the corresponding Chebyshev polynomial will converge. Using the Newton formula, this 
result of convergence is true only on the theoretical level. It is shown that the algorithm 
which computes the divided differences is numerically stable only if: 1.) the interpolating 
points are arranged in a certain order, 2.) the size of the interval is 4. 
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1. INTRODUCTION 


Let f(x) be a real function defined on the interval [a, 6] and {x,-}£L 0 be a set of N + 1 
points in [a, 6] then the general formulation of Newton interpolating polynomial of degree 
N is 

N 

Pn{x) = J2 a k R k (x) (1.1) 

Jt=0 

where a k are the divided differences 


a k = f[x o, — , as*] (1-2) 

and 

R 0 {x) = 1 (1.3) 

R k (x) = - x { ). (i.4) 

If En(x) is the error at point xe[a, 6] then we have the following theorem [CodB82], 


Theorem 1.1: Let f(x) be a real valued function defined on [a, 6] and N + 1 times 
differentiable on [a, b\. If Pn[x) is the polynomial of degree < N which interpolates f(x) 
at IV + 1 distinct points x 0 , • • • ,xn in [a, 6], then for all xe[a,b], there exists f = £(x)e(a, b ) 
such that 

«»(*) = m - jvw = (i-5) 

It is well-known that if xq, • • • ,x^ are equally spaced points then max\RN + i(x)\ 
increases as x moves towards the ends of the interval and divergence can occur. A well- 
known remedy for this phenomenon is to choose as interpolating points the zeroes of the 
corresponding scaled and translated Chebyshev polynomial 

Xi = \ (6 ~ a)cos + b + a i = 0,-",N. (1.6) 

Using (1.6) results in a uniform distribution of the error and convergence is achieved. 

This result of convergence is true only in theory. Practically, an interpolating polyno- 
mial based on (1.6) will not converge to f(x) because of the finite accuracy of the computer. 
The numerical instability can be traced to two sources: 
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1. ) The algorithm which computes the divided differences is very vulnerable to roundoff 

errors and overflow. The super sensitivity to roundoff errors is explained by the fact that 
for N large, the first k points ( k N) are very close to each other (mm|x i+1 — x,\ ~ 1/N 2 ). 

But even if roundoff errors were eliminated, we still would face overflow because the first 
k points ( k <C N) are concentrated on one side of the interval. This distribution will lead 
to nonuniformity of R k (x) and E k {x). It is this nonuniformity which will cause overflow at 
intermediate stages of the interpolation process. 

2. ) Eliminating the first source of numerical instability (by taking the points in a 
different order) results in “almost” uniformity of R k {x ) k — 1, .., N, but still, R k (x) satisfies 
(see Section 3) 

\R k {x)\~2~ k . (1.7) 

Observing (1.1) and (1.7) it is obvious that we will face overflow (for k large enough) while 
computing the a' k s. 

We are going to approach this phenomenon in the more general context of interpolation 
in the complex plane. Background material is given in Section 2. Based on the theory we 
show that by: 

(1) arranging the interpolating points in a certain order, 

(2) making a simple change of variables such that the interval is of size 4, 
we get 

|i4(x)|~l 1 <k<N. (1.8) 

Thus these two modifications result in a stable Newton interpolation process. Algorithm 
2.1 presented in Section 2 generates N interpolating points such that (1.8) is satisfied. 
In practice we would like to be able to add points, if necessary, without restarting the 
interpolation process. It is well-known that the Newton algorithm potentially has this 
feature; adding another point results in computing only one additional term. (In contrast 
to the Lagrange process where one has to start the calculations all over again.) But by 
using Chebyshev zeroes, (1.6), we would not capitalize on this property since increasing 
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N results in changing all the previous points. To this end, we construct Algorithm 2.2 
which enables us to add interpolating points such that (1.8) is asymptotically satisfied. 
Observe from (1.1) and (1.8) that the a' k s behave asymptotically like the interpolation 
error. Thus one can use the algorithm which computes the divided differences (related to 
the points generated by Algorithm 2.2) as a numerical device for estimating the degree 
of the polynomial needed to achieve a given accuracy. This technique can be very useful 
when the mathematical expression of the error (1.5) is difficult to analyze. 

The need for high degree interpolation is confronted in approximating a finite operator 
which can be presented as f(A) where A is another finite operator. (The case f(z) = exp(z) 
is a popular example [MoVa78].) Approximating f(A) can be reduced to a problem of 
approximating f(z) where z belongs to a domain D in the complex plane which includes all 
the eigenvalues of A. A possible approach is to expand the function as a sum of orthogonal 
polynomials. In [Tale86], [Tale85] we have used it for the function f(z ) = exp(tz ) where 
the domains were [-iR,iR], [-R,0] respectively. The algorithm which results make use of 
the three term recurrence relation satisfied by the polynomials (scaled and translated 
Chebyshev). For more complicated domains, the related polynomial of degree k satisfies 
a k term recurrence relation and therefore the expansion approach is not suitable. An 
alternative way is to use interpolation. A brief description of the method is given in 
Section 4. 

A particular and very important case of polynomial approximation of / (A) is an itera- 
tive solution to a linear system Ax = b. Finding optimal parameters a* for solving Ax = b 
by the Richardson algorithm 

x k+1 = x k - cc k {Ax k - b) (1.9) 

can be achieved by considering polynomial interpolation to the function f(z) = 1 /z 
[Tale87]. When D is on one side of the real line, this method is widely treated in the 
numerical analysis literature (Chebyshev acceleration [HaYo81]). It is also known that the 
fact of having to decide on the number of iterations before starting the algorithm reduces 


3 



its efficiency. A Richardson process which uses points generated by Algorithm 2.2 is free of 
this disadvantage. (For a more elaborate discussion and numerical examples see [Tale87].) 
We conclude the paper in Section 5 by giving some numerical results. 


2. INTERPOLATION IN THE COMPLEX PLANE 


Let D be a bounded continuum in G such that the complement of D is simply connected 
in the extended plane and contains the point at infinity. Considering interpolation in the 
domain D, one is faced with the problem - which are the “good” interpolating points? 
The solution to this question is based on the following: 

Let <f>(z ) be a conformal mapping which maps the complement of D to the complement 
of a disc of radius p such that 


Hm M = !. 


( 2 . 1 ) 


z— ► OO £ 

p is the logarithmic capacity of D [SmLe68]. (Having a domain D,<f> and p are defined 
uniquely.) Define ip(oj) to be the inverse of (f>(z ). Then we have [Wals56]: 


Definition 2.1: Let Tjj be the image under rp of the circle |iu| = R (R > p) and Ir 
be the closed Jordan region whose boundary is T^. If f(z) is single valued and analytic on 
Ir then the sequence of polynomials P m (z) is said to converge to f(z) on D maximally if 

\f(z)-P m (z)\<C(p/R) m zeD (2.2) 

where C depends on p/R but not on m or z. 

Definition 2.2: The set of interpolating points Zj = ip(wj) is said to be uniformly 
distributed on Td (the boundary of D ) if wj are equally distributed on the circle \w\—p. 

This set (known also as Fejer points) is one possible set of “good” points. The polyno- 
mial which interpolates f(z ) at these points satisfies (2.2) [Wals56]. 

Another possibility is to interpolate at the zeroes of the corresponding Faber polyno- 
mial. Let D and <p(z) be as defined previously. We have [Mark77] 

[<j>{z)] m = C m z m + • • - + C 0 + C-i/z + C. 2 /z 2 + ■■■. (2.3) 
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The Faber polynomial of degree m related to D is the polynomial part of (2.3) 

F m (z) = C m z m + --- + C 0 . (2.4) 

Interpolating at the zeroes of F m (z) satisfies (2.2) [Mark77]. When 

D = {z\ — 1 < Rez < 1 I m z = 0}, (2.5) 

then 

F m {z) = T m (z) = cos(mcos _1 ( 2 )), (2.6) 

and the m zeroes are the Chebyshev points (1.6). 

For a general domain in the complex plane finding the zeroes of F m (z) for large m can 
be troublesome. Thus, it is preferable to use Fejer points since only knowledge of ijj{w ) is 
required. 

Assume now that Zj, 1 < j < m, are uniformally distributed points. Then [Wals56]: 

\Rm(z)\ ~P m . (2.7) 

From (2.7) it is clear that in order to satisfy (1.8) we need: 

1. p = 1 

2. Every subset of interpolating points 2 ,-, • • • , Zk, 1 < k < N, has to be uniformly 
(or “almost” uniformly, see Algorithm 2.1) distributed. 

Hence, by making the change of variables 

z = z/ p (2.9) 

and arranging the interpolating points such that the second requirement is satisfied we 
eliminate the numerical instability mentioned in Section 1. The following algorithm is 
designed for this purpose. 

Algorithm 2.1: [Generates m uniformly distributed points (m even)] 

6i = 0 j = t = 1 86 = 2tt /m. (2.13) 
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Find the largest k such that k is power of 2 and k < mj 2 


6k = k x 86 ( 2 -14) 

1 For i — 1 until i do: 

r] = 6i + 0A; 

*/(*? > ?r) go to 1 

i = 3 + 1 (2.15) 

Oj = t) 

end do 

t = j 
6k = 6k/2 
if {Ok < 60) stop 
go to 1 

Algorithm 2.1 generates y arguments of points on the upper part of the unit circle (includes 
1 but not -1). Thus 

u>i = 1 Zi = ip(w i) (2.16a) 

w 2 = -1 z 2 = 4 >{w 2 ) (2.166) 

771 

w 2 j-i = exp(i0j) z 2 j-i = r/j(w 2 j-i) 2 < j < — (2.16c) 

77 % 

w 2 j = —w 2 j-i z 2j - ip(wij ) 2 <j<~. (2.16d) 

Using Algorithm 2.1 results in uniform distribution of Zi, • • • , z m while z lf • • • , zt(k < m) 
are “almost” uniformly distributed. 

As mentioned in the introduction, we would like to h?,ve an interpolating process which 
will allow us to add points if desired. To this end, we present the next algorithm. It 
generates an infinite set of interpolating points. Using this algorithm we do not have 
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decide on the degree of the polynomial ahead of time. The set is asymptotically uniformally 
distributed and therefore the interpolating polynomial satisfies (2.2). 

Algorithm 2.2: 

60 = 7T 

k = 1 
z\ = V’(l) 

1 For i = 1 until k or until satisfied do 

Ok+i = 0i + 60 
Zk+i = ip[exp(iO k+i )) 

end do 

86 - 60/2 
k = 2k 
go to 1 

The set of points generated by Algorithm 2.2 is uniformally distributed when the num- 
ber points is a power of two and “almost” uniformally distributed otherwise. 

Remark : In order to implement the algorithms presented in this section, we need the 
conformal mapping function. For some domains, we do have an analytic expression of 
this function. In more complicated situations, one has to resort to numerical techniques. 
When D is a polygon, ip(w) is a Schwartz-Cristoffel transformation. Numerical routines 
for this case, written by L. N. Trefethen (based on [Tref80]), are available through the 
Netlib facilities. (See also [Tale87] for a description of how to implement the routines.) 
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3. INTERPOLATION ON THE REAL LINE 


Let 

D = [-b,b}. (3.1) 

The conformal mapping xp(w) which maps the complement of the unit disc on the comple- 
ment of D is [Mark77] 


xj)( w) = -(to + xv x ). 

Using (2.1) we get that 

p = 6 / 2 . 

Therefore, one should make the following change of variables 

.. 2x 


x = 


in order to get 


Similarly, if 


then 


£ = 1 - 2 , 2 ] 

P= I- 

D = [a, 6] 


x = 


4.x 

6 — a 


(3.2) 

(3.3) 

(3.4) 

(3.5) 

(3.6) 

(3.7) 

(3.8) 


Thus, without loss of generality, we assume that 


D = j—2, 2]. 


(3.9) 


Using (3.2), we get that Fejer points are 

Xj = xv j + xv J 1 (3.10) 

where Wj can be generated by the algorithms given in Section 2. Observe that the Xj 's 
are double interpolation points (except for 2 and -2). 
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Faber polynomials which correspond to D are scaled Chebyshev polynomials [Mark77] 

T m (z) = cos(m(cos _1 (|))) (3.11) 

whose zeroes are 

(2j — lW . . . 

yj = 2cos - — j = 1, • • • ,m. (3.12) 

2m 

Using Algorithm 2.1, we can get the zeroes of T m (x), arranged in a stable order, as follows: 

60. 

Xj = 2 cos{Qj + — ) j = 1, • • • , m. (3.13) 

A popular set of interpolating points on the real line is the extremas of Chebyshev 
polynomial of degree N. This set is not exactly Fejer points but can be shown to satisfy 
(2.2). Using the algorithms given in Section 2 we have 

X\ = 2; x 2 = —2 (3.14a) 

Xj = 2 cosOj j > 3 (3.146) 

where Oj are generated by Algorithm 2.1. or 

Xi = 2 x 2 = —2 (3.15a) 

Xj = 2cos$2j+i j > 1 (3.156) 

where 0 2 j+i are generated by Algorithm 2.2. 

4. POLYNOMIAL APPROXIMATION OF A FUNCTION OF A MATRIX 

Let A be an N x N matrix and f(z) a function analytic in a domain D in the complex 
plane which includes all the eigenvalues of A. Let w be the vector which results from 
operating with f{A) on a vector v 

w = [fi A )\ v - (4-1) 

Getting an approximation of w is a problem frequently confronted in applied mathematics. 
An elaborate description of this topic is given in [Tale87]. It is shown there that one 
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can get an “almost” optimal polynomial algorithm by using the polynomial P m {z) which 
interpolates f(z ) at points uniformally distributed on the boundary of D. (Without loss 
of generality we assume that p(D) = 1. If not, we define A — - A and consider f(z) = 
f(pz) = f(z)). Let vj be an approximation of w 

w = [P m (A)]v. (4.2) 

When A is normal, the error vector satisfies 


\ W ~ M\ = ll(/U) - Pm{A))v\\ < | f{z) - P m {z)\ |M|. (4.3) 


Hence 


E = 


|u; — w\ 

|u| 


< | f{z)-P m {z)\. 


(4.4) 


Using (2.2) we get that E is asymptotically bounded by (|>) m . Thus one can get a prediction 
of the degree of the approximating polynomial. (When f(z) is an entire function ( R = oo) 
one can use the idea of computing the corresponding divided differences in order to get an 
estimation of m.) In this case, we will use Algorithm 2.1 as a generator of the interpolating 
points Zi. Once we have z, we compute the divided differences a,-. Having these two sets of 
numbers, we can write the following algorithm: 


Algorithm 4.1: 

u = v 

W = CliU 

For * = 2, • • • , m do 

u — (A — z,_i/)u 

W — W + diU 

end do 

(4.5) 


When A is diagonalizable but not normal we have the following 


E < imi IIT-'ll \pz)-p m (z)\ 
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where T is the matrix which diagonalize A. Since ||T|| ||T -1 || is unknown we usually 
would not have an estimation of m before starting the algorithm. In this case we can use 
Algorithm 2.2 as a generator of points. The criterion for stopping will be the norm of the 
relative residual ||a,u||/||tl/||. (Because of the asymptotic behavior, the decision should be 
made checking by a few residuals.) 

Algorithm 4.2: 

u = v 
w = aiu 

For i = 2, • • •, until satisfied do 

u = (A — Zi_il)u 

W = W + dill 

check for convergence 

When A is a real matrix and f(z) is real for z real, the algorithms have to be modified in 
order to eliminate complex arithmetic. (Observe that in this case the domain is symmetric 
around the real axis.) To this end, we first arrange the points such that each point with 
nonzero imaginary part will be followed by its conjugate. Thus, we change (2.16d) to read 
Wtj = uT 2 j _i. The modified version of Algorithm 4.1 is: 

Algorithm 4.3: 

w = [ ail + a 2 (A — zil)]v 
r = (A — ziI)(A — z 2 I)v 

For * = 1, • • • , do 

f= (A- z? i+1 I)r 
w = w + a^r + a£ +1 r 
r = (A - z? i+1 I)f + (4 + i) 2 r 

end do 
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(The superscripts R, I stand for the real and imaginary parts respectively.) 


Using Algorithm 2.2, we will order the points as follows: z 2 , 23, z 3 , • • •, where z 2 »+i 

are generated by the algorithm. The modified version of Algorithm 4.2 will read: 

Algorithm 4.4: 

w = [axJ + a 2 (A — Zil)\v 
r = (A — ziI)(A — z 2 I)v 

For i = 1, • • •, until satisfied do 

r = (A — z 2i+1 I)r 
w = w + a*r + a£ +1 r 
r = {A- z* +1 I)r + (4 + i) 2 r 
check for convergence 


An iterative solution to a general linear system 

Ax = b (4.7) 

is a particular case of approximating f(A) where f(z) = K Applying Algorithms 4.1 - 
4.4 for this problem can be shown to be equivalent to a Richardson type process [Tale87]. 
Having an initial approximation x 1 , Algorithm 4.1 will result in 

Algorithm 4.5: 

For j = 1 until m do 

x j+1 = x i — aj(Ax* — b ) 

end do 


Where the relaxation parameters are 


a i = — • 


(4.8) 
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Thus, for this particular function we do not have to compute the divided differences. 
Therefore, the change of variables (2.9) is not significant here. But still rearranging the 
points is important. Otherwise, we can face numerical instability due to the nonuniformity 
of the error at intermediate stages. This phenomenon is treated also in [AnGo72], [LeFi76], 
and [FiRe87]. The algorithm (based on bit-reversed binary representation of the iteration 
count) proposed in [FiRe87] results in the same set of parameters generated by Algorithm 
2.2. When A is real, Algorithm 4.4 will read 

Algorithm 4.6: 

x 2 = x 1 — a.\(Ax l — b ) 
x 3 = x 2 — a 2 (Ax 2 — b) 

for j = 2 until satisfied do 

R’ = Ax 2 ’- 1 - b 

x 2 ’ +1 = X 2 ’- 1 - [2 af - \ocj\ 2 A\R’ 
check the residuals 

end do 

5. NUMERICAL RESULTS 


All the numerical experiments reported in this section were performed on an IBM 3270. 
Let P m (x) be the interpolating polynomial. The maximal error 


E= max |/(x) 
*« [-i.il 

Pm(x) | 

(5.1) 

was approximated by 



E ap = max \f(xi) 

Jm(®«)| 

(5.2) 

where £,• are checkpoints 



x { = — 1 + 2(* — 1)/19 

1 < * < 20. 

(5.3) 
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In the first set of experiments, the interpolated function is 


f( x ) = 


1 

0.005 + x 2 


~ 1 < x < 1 . 


By the following change of variables 


(5.4) 


y = xb 


(5.5) 


we get 

fix) = giy) — — — — — b < y < b. (5.6) 

v ) 0.005+ [I) 2 ~ ~ K ’ 

The interpolating points are uniformally distributed on [—6,6] (as defined in Section 2) 

but taken in two different orderings: yj are arranged in the standard order 

y) = 6cos(7rj/m) 0 < j < m (5.7) 

and yj are generated by Algorithm (2.1). In the first experiment we have used single 
precision. The polynomial is of degree 80. Selected divided differences and the numerical 
maximal error are presented in Table 1. 


Table 1. (m = 80) 


k 

b- 

=1 

b= 

-2 

b= 

=3 


ak{y l ) 

Gfc(y 2 ) 

^(y 1 ) 

ak[y 2 ) 

My 1 ) 

My*) 

Hul 

-5.44+09 

-4.75+03 

-1.06+07 

-9.28+00 

-2.08+04 

-1.81-02 

1 

-5.31+17 

2.045+06 

-1.02+12 

3.89+00 

-1.93+06 

7.44-06 


-1.75+21 

-8.94+08 

-3.26+12 

-1.66+00 

-6.07+03 

-3.10-09 

40 

-1.46+22 

3.45+11 

-2.65 +10 

6.36-01 

-4.83-02 

1.16-12 

50 

-2.54+21 

-2.14+14 

-4.50+06 

-3.80-01 

-8.00-09 

-6.74-16 

60 

-4.28+21 

9.27+16 

-7.42+03 

1.61-01 

-1.29-14 

2.80-19 

70 

-1.73+22 

-3.43+19 

-2.94+01 

-5.82-02 

-4.98-20 

-9.86-23 

80 

-8.79+16 

2.90+18 

-1.46-07 

4.80-02 

-2.41-31 

7.94-30 

E ap 

4.68+23 

3.48-01 

4.68+23 

3.48-01 

— 

3.48-01 


Table 1 verifies the fact that yj are useless for high degree Newton interpolation. Since 
no overflow (underflow) has occurred in computing the divided differences, the error related 
to y ? is the same in all the three cases. However, observe the “nice” behavior of the 
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coefficients which correspond to b = 2. The last ones are close to the interpolation error, 
as predicted by the theory. 

Next we took m = 160. The results are presented in Table 2. 


Table 2. (m = 160) 


k 

b= 

=1 

b= 

=2 

b= 

=3 


My 1 ) 

My 2 ) 

My 1 ) 

ajt(y 2 ) 

My 1 ) 

My 2 ) 

20 

3.91+29 

2.42+06 

7.46+23 

4.62+00 

1.42+18 

8.81-06 

40 

— 

5.44+11 

6.49+33 

9.90-01 

1.18+22 

1.78-12 

60 

— 

1.36+17 

9.31+34 

2.35-01 

1.65+17 

4.08-19 

80 

— 

2.34+17 

7.43+30 

3.88-02 

1.23+07 

6.43-26 

100 

— 

7.24+27 

2.10+23 

1.14-02 

3.31-07 

1.80-32 

120 

— 

1.67+33 

3.28+13 

2.51-03 

4.93-23 

3.77-39 

140 

— 

— 

1.01+02 

2.45-04 

1.45-40 

0.00+00 

160 

— 

— 

-5.35-08 

-2.60-05 

0.00+00 

0.00+00 

E ap 

— 

— 

— 

7.29-04 

— 

— 


Now, overflow or underflow is taking over in all the cases except for y? and 6 = 2. As 
in the previous example, the last coefficients are close to the error itself. 

Next we interpolated the function 


f{x) = cos(2000x) - 1 < x < 1. 


(5.8) 


In this case, one needs a super high degree polynomial. We took m = 2100 and used 
double precision. The divided differences are presented in Table 3. 


Table 3. (m = 2100) 


k 

1 M 

k 

IM 

k 

M 

k 

IM 

1 

3.67-01 

801 

2.36-02 

1601 

1.71-02 

2031 

1.98-05 

201 

8.52+00 

1001 

3.80-03 

1801 

4.47-02 

2051 

1.12-06 

401 

3.78-01 

1201 

1.51-02 

2001 

5.13-02 

2071 

5.77-08 

601 

5.69-03 

1401 

1.22-02 

2011 

2.78-03 

2091 

2.18-10 


The numerical error is 


E ap = 5.89 - 09. (5.9) 

Observe that the asymptotic decay of the coefficients starts when k crosses 2000. One 
encounters a similar behavior of coefficients while expanding (5.8) in Chebyshev polyno- 
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mials. 


OO 

cos(2000x) = Yl b k T k {x). - 1 < x < 1 (5.10) 

k=0 

It can be shown [AbSt72] that 


bk { 2Jjt(2000) k even ( 5 ' U ) 

where J k are Bessel functions of order k. By asymptotic analysis of Bessel functions we 
know that J k {R) goes to zero exponentially fast only when k > R. This similarity between 
a k and b k could be anticipated since T k [x) satisfies 

1. |T*(*)|<1 |*| <1 

2. Its zeroes are uniformally distributed. 

Where R k (x) satisfies 

1. |Jfc(*)| ~ 1 |*| <2 

2. Its zeroes are “almost” uniformally distributed. 


Polynomial approximation of a function of a matrix 
Let us consider the following parabolic P.D.E 


U t = U xx + U vy |*| < 1, |y| < 1 

4 

U(x, 0) = 5>n(***)«»(jrlfcy) 

Jb=l 


U = 0 on the boundaries of the unit square. 
Using the following space discretization 


(5.12a) 

(5.125) 


Uii = U( Xi 


(5.13 a) 

Xi = — 1 + iAx 

1 <i<N 

(5.135) 

Vj = -1+jAy 

1 <j<N 

(5.13c) 

> 

H 

II 

> 

II 

2 

(5.13d) 

N + l 
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where N is the number of grid points in each direction, we approximate the space deriva- 
tives by the standard second order central differencing 


U xx ~ {U i+1>j - 2 Uij + U { _ 1,,0/Ax 2 (5.14a) 

Uyy ~ (Ui, j+ 1 - 2C7,y + Uij-i)/ Ay 2 . (5.146) 

Thus we get the following semidiscrete representation of (5.12) 

(t^iv)t = GnUn (5.15a) 

= 17° (5.156) 


where Z7jv is a vector of dimension iV 2 , is a N 2 x iV 2 matrix which results from the 
space discretization and 

Cl?. = y y , 0) 1 < t, j < N. (5.16) 

The formal solution of (5.15) is 

U N = exp{tG N )U° N (5.17) 

Now, by using polynomial approximation of exp(tG ^ ) we get 

U k N = P k {tG N )U° N (5.18) 

where P k (tz) approximates exp(tz) in the domain which includes all the eigenvalues of Gjy. 
The numerical solution was computed at t = 0.1. A simple Fourier analysis shows that 
D is on the negative real line. Therefore we took (3.126) as interpolating points and used 
Algorithm 4.2 to get U^. The exact solution of (5.12) is 

4 

U(x,y,t) = T, exp{—[k-K) 2 t)sin[k'nx)sin[k'ny). (5.19) 

k=l 

The following tables present the dependence of three factors on k, the degree of the inter- 
polating polynomial. 
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|<zjt| — the absolute value of the maximal divided difference, 
fjfc/ujfc— the relative residual which represents the time accuracy. 
e k/ u k~ the relative error which represents the space + time accuracy. 
For N = 8 D is 


D = {x\ - 160 < x < -4} (5.20) 

and the results are given in Table 4. 


Table 4. ( N = 8) 


k 

M 

rk/u k 

e k/ u k 

k 

l°*l 

rk/u k 

e*/ujfc 

2 

2.45-01 

7.79+00 

7.01+00 

12 

2.73-03 

1.85-02 

8.30-02 

4 

1.79-01 

5.42+00 

3.27+00 

14 

2.25-05 

1.92-03 

8.23-02 

6 

8.29-02 

2.32+00 

1.74+00 

16 

1.11-06 

1.48-04 

8.23-02 

8 

1.31-02 

5.81-01 

2.54-01 

18 

6.15-08 

1.77-06 

8.23-02 

10 

2.40-03 

9.18-02 

1.06-01 

20 

2.97-09 

1.94-07 

8.23-02 


For N — 16 D is 


D = {x | - 640 < x < -4} (5.21) 


and the results are given in Table 5. 


Table 5. (JV = 16) 


k 

O-k 

rk/uk 

ek/u k 

k 

a.k 

rk/u k 

Ck/uk 

3 

1.67-01 

2.44-01 

9.19+00 

18 

9.50-04 

3.28-02 

9.76-02 

6 

2.88-02 

1.67-01 

2.64+00 

21 

6.32-05 

1.80-03 

2.28-02 

9 

2.29-02 

5.19-01 

3.59+00 

24 

1.64-05 

3.05-03 

2.29-02 

12 

4.08-02 

1.35+00 

2.36-01 

27 

5.92-06 

3.71-04 

2.26-02 

15 

6.50-03 

3.05-01 

2.46-01 

30 

2.56-07 

4.96-05 

2.26-02 


Observing Tables 4 and 5 we see that for N = 8 we need 14 iterations to recover the 
space accuracy and for N — 16 we need 27 iterations. This result compares favorably with 
a standard explicit algorithm like Forward Euler or Runge-Kutta where the number of 
iterations(time steps) is multiplied by 4 (due to stability limitations) whenever we increase 
the space resolution by two. In [Tale85] we describe an algorithm which is spectral in time. 
It is shown there that for a parabolic problem, the time resolution parameter depends 
almost linearly on the space resolution parameter. Since the algorithm used in the present 
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paper can be considered also as spectral in time, it explains this improvement in the number 
of matrix vector multiplications needed to march the solution to a given time level. 
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